home *** CD-ROM | disk | FTP | other *** search
/ Meeting Pearls 1 / Meeting Pearls Vol 1 (1994).iso / installed_progs / gfx / lise2.1 / lise / src / mspc.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-03-31  |  10.8 KB  |  383 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include <spec.h>
  4.  
  5. int umax,ulen,ubeg,uend;
  6. float *spc,*err,*tim,*spcx,*errx,*timx;
  7.  
  8. char *comment,*uname;
  9.  
  10.  
  11. help()
  12. {
  13. printf("modify spectrum\n");
  14. printf("mspc sp{[i:j]} n1 n2 ... {-sin} {-cos} {-exp} {-ln} {-sqr} {-sqrt} {-abs}\n");
  15. printf("                         {-s n} {-d n} {-p n} {-q n} \n");
  16. printf("the contents of sp is modified starting from i, ending with j. \n");
  17. printf("contrary to the normal behaviour, the WHOLE spectrum (and corresponding\n");
  18. printf("error array) is stored back \n");
  19. printf("If n1 n2 ... is specified, then sp[i:j] is set to the constant values \n");
  20. printf("with the last one repeated up to j. \n");
  21. printf("other operations: \n\n");
  22.  
  23. printf("   -f   use as Filter (do not write back) \n\n");
  24.  
  25. printf("   -c text    change comment only. \n");
  26. printf("   -tica n.m  change time calibration for this file \n");
  27. printf("   -t    modify timespectrum only. \n");
  28. printf("   -e    modify errorspectrum only. \n\n");
  29.  
  30. printf("   -sin  for all n: spc(n) = sin(spc(n)) \n");
  31. printf("   -cos  for all n: spc(n) = cos(spc(n)) \n");
  32. printf("   -exp  for all n: spc(n) = exp(spc(n)) \n");
  33. printf("   -ln   for all n: spc(n) = log(spc(n)) \n");
  34. printf("   -sqr  for all n: spc(n) = spc(n)*spc(n) \n");
  35. printf("   -sqrt for all n: spc(n) = sqrt(spc(n)) \n");
  36. printf("   -abs  for all n: spc(n) = fabs(spc(n)) \n");
  37. printf("   -int  for all n: spc(n) = floor(spc(n)) \n");
  38. printf("   -lin  for all n: spc(n) = n \n");
  39. printf("   -mir  for all n: spc(n) = spc(max-n), modifies error array too \n\n");
  40.  
  41. printf("   -add x  for all n: spc(n) = spc(n) + x\n");
  42. printf("   -sub x  for all n: spc(n) = spc(n) - x\n");
  43. printf("   -mul x  for all n: spc(n) = spc(n) * x\n");
  44. printf("   -div x  for all n: spc(n) = spc(n) / x\n\n");
  45.  
  46. printf("   -sr  x for all n: spc(n+x) = spc(n)  (shift right) \n");
  47. printf("   -srt x for all n: spc(n+x) = spc(n)  (shift right, including time) \n");
  48. printf("   -sl  x for all n: spc(n) = spc(n+x)  (shift left) \n");
  49. printf("   -slt x for all n: spc(n) = spc(n+x)  (shift left, including time) \n");
  50. printf("   -com x compress spectra, modifying error and timearrays too\n");
  51. printf("   -bel x smoothly cut spectrum after x, modifies error array too\n\n");
  52.  
  53. printf("instead of numbers, other spectra arrays may be specified. \n");
  54. printf("these spectra arrays are taken starting from the first element.\n");
  55. printf("ONLY ONE OPERATION CAN BE PERFORMED AT A TIME ! \n");
  56. return(0);
  57. }
  58.  
  59. main(argc,argv)
  60. int argc;
  61. char *argv[];
  62. {
  63. int   m,n,i;
  64. int   max;
  65. char  *s,*z;
  66. FILE  *stream,*fopen();
  67.  
  68.  
  69.    if(argc<3) _help0(); /* at least two parameters should be specified */
  70.  
  71.    z = (char *) malloc(80);
  72.    s = (char *) malloc(80);
  73.    comment = (char *) malloc(80);
  74.    uname = (char *) malloc(80);
  75.  
  76.    checkopt(argc,argv,"-dummy",z);
  77.  
  78.    spc= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  79.    err= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  80.    tim= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  81.    spcx= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  82.    errx= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  83.    timx= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  84.    if(timx==NULL) {
  85.       printf("sorry, not enough memory\n");
  86.       exit(-1);
  87.    }
  88.  
  89.    if(checkopt(argc,argv,"-lin",z)) { /* this should be possible, without
  90.                                          having a spectra array yet */
  91.       for(n=1;n<=argc;n++) {
  92.          i=instr("[",argv[n]);
  93.          if(i>0) {
  94.             midstr(z,argv[n],0,i-1);
  95.          }
  96.       }
  97.       strcpy(s,z);
  98.       stream=fopen(z,"r");
  99.       if(stream==NULL) {
  100.          strcat(z,".spc");
  101.          stream=fopen(z,"r");
  102.          if(stream==NULL) {            /* create a minimum spectrum */
  103.             strcpy(z,"himmelarschundzwirn");
  104.             writespec(s,spc,err,2,0,z);
  105.          }
  106.       }
  107.    }
  108.  
  109.    
  110.    max=readspec(argv[1],spc,err,tim,comment);
  111.    umax=_max; ubeg=_sbeg; uend=_send;
  112.    strcpy(uname,_spc_0nam);
  113.    if(checkopt(argc,argv,"-o",z)) strcpy(uname,z); /* use another output */
  114.    for(n=0;n<=umax;n++) {
  115.       spc[n]=_uspc[n];
  116.       err[n]=_uerr[n];
  117.       tim[n]=_utim[n];
  118.    }
  119.  
  120.    if(checkopt(argc,argv,"-f",z)) strcpy(uname,""); /* use as filter */
  121.    if(checkopt(argc,argv,"-tica",z)) { /* modify time cal. */
  122.       _tica = atosf(z);
  123.       for(n=0;n<=umax;n++) tim[n]=_tica * n;
  124.    }
  125.    if(checkopt(argc,argv,"-t",z)) { /* modify time spectrum only */
  126.       strcat(uname,".tim");
  127.       for(n=0;n<=umax;n++) spc[n]=tim[n];
  128.    }
  129.    if(checkopt(argc,argv,"-e",z)) { /* modify error spectrum only */
  130.       strcat(uname,".err");
  131.       for(n=0;n<=umax;n++) spc[n]=err[n];
  132.    }
  133.  
  134.    doit(argc,argv);
  135.  
  136.    writespec(uname,spc,err,umax,2,comment);
  137.    if(instr(".tim",uname) < 0) {
  138.       strcat(uname,".tim");
  139.       writespec(uname,tim,err,umax,2,comment);
  140.    }
  141.  
  142.    free(spc); free(err); free(tim); 
  143.    free(spcx); free(errx); free(timx);
  144.    free(uname); free(comment); free(z); free(s);
  145.    exit(0);
  146. }
  147.  
  148. doit(argc,argv)
  149. int argc;
  150. char *argv[];
  151. {
  152. float x,y;
  153. int   m,n,i,j,k,max;
  154. char  c,z[80];
  155.  
  156.    for(n=1;n<=argc;n++) {               /* check for comment modification */
  157.       if(strcmp(argv[n],"-c")==0) {    /* modify comment */
  158.          strcpy(comment,"");
  159.          for(i=n+1;i<=argc;i++) {
  160.             strcat(comment,argv[i]);
  161.             strcat(comment," ");
  162.          }
  163.      return(0);
  164.       }
  165.    }
  166.  
  167.    max=0;                              /* search for number or spectra name */
  168.    for(n=2;n<argc;n++) {
  169.       c=argv[n][0];
  170.       if(c=='-') {
  171.          c=argv[n][1];
  172.          if(c>'9') break;              /* options reached, no number or file name */
  173.          if(strlen(argv[n])==2) break;
  174.          spcx[max++]=atof(argv[n]);     /* assume negative number */
  175.       }
  176.       if((c>='0') && (c<='9')) {       /* catch number(s) */
  177.          spcx[max++]=atosf(argv[n]);     /* assume number */
  178.      continue;
  179.       }
  180.       max=readthis(argv[n],spcx,errx,timx); /* possible spectra array */
  181.    }
  182.  
  183.    if(max>0) {                         /* set to constant */
  184.       for(i=0;i<max;i++) spc[ubeg++]=spcx[i];
  185.       while(ubeg <= uend) spc[ubeg++]=spcx[i-1];
  186.       return(0);
  187.    }
  188.  
  189. /* ----------------- binary functions ------------------ */
  190.       
  191.    if(checkopt(argc,argv,"-add",z)) {    /* sum */
  192.       max=getnums(z,spcx,errx,timx);
  193.       for(i=0;i<max;i++) {spc[ubeg] = spc[ubeg] + spcx[i]; ubeg++ ;}
  194.       while(ubeg <= uend) {spc[ubeg] = spc[ubeg] + spcx[i-1]; ubeg++ ;}
  195.       return(0);
  196.    }
  197.    if(checkopt(argc,argv,"-sub",z)) {    /* difference */
  198.       max=getnums(z,spcx,errx,timx);
  199.       for(i=0;i<max;i++) {spc[ubeg] = spc[ubeg] - spcx[i]; ubeg++ ;}
  200.       while(ubeg <= uend) {spc[ubeg] = spc[ubeg] - spcx[i-1]; ubeg++ ;}
  201.       return(0);
  202.    }
  203.    if(checkopt(argc,argv,"-mul",z)) {    /* product */
  204.       max=getnums(z,spcx,errx,timx);
  205.       for(i=0;i<max;i++) {spc[ubeg] = spc[ubeg] * spcx[i]; ubeg++ ;}
  206.       while(ubeg <= uend) {spc[ubeg] = spc[ubeg] * spcx[i-1]; ubeg++ ;}
  207.       return(0);
  208.    }
  209.    if(checkopt(argc,argv,"-div",z)) {    /* quotient */
  210.       max=getnums(z,spcx,errx,timx);
  211.       for(i=0;i<max;i++) {
  212.          if(spcx[i]!=0.0) spc[ubeg] = spc[ubeg] / spcx[i];
  213.          ubeg++ ;
  214.       }
  215.       while(ubeg <= uend) {
  216.          if(spcx[i-1]!=0.0) spc[ubeg] = spc[ubeg] / spcx[i-1];
  217.          ubeg++ ;
  218.       }
  219.       return(0);
  220.    }
  221.  
  222. /* -------------------- shifting and compression ---------------- */
  223.  
  224.    if(checkopt(argc,argv,"-sr",z)) {    /* shift right */
  225.       m=atoi(z);
  226.       for(i=uend-m;i>=ubeg;i--) {
  227.          spc[i+m] = spc[i];
  228.          err[i+m] = err[i];
  229.       }
  230.       for(i=ubeg;i<m;i++) {
  231.          spc[i] = 0.0;
  232.          err[i] = 0.0;
  233.       }
  234.       return(0);
  235.    }
  236.  
  237.    if(checkopt(argc,argv,"-srt",z)) {    /* shift right with time*/
  238.       m=atoi(z);
  239.       for(i=uend-m;i>=uend;i--) {
  240.          spc[i+m] = spc[i];
  241.          err[i+m] = err[i];
  242.          tim[i+m] = tim[i];
  243.       }
  244.       for(i=ubeg;i<m;i++) {
  245.          spc[i] = 0.0;
  246.          err[i] = 0.0;
  247.          tim[i] = 0.0;
  248.       }
  249.       return(0);
  250.    }
  251.    if(checkopt(argc,argv,"-sl",z)) {    /* shift left */
  252.       m=atoi(z);
  253.       for(i=ubeg;i<=uend;i++) {
  254.          spc[i] = spc[i+m];
  255.          err[i] = err[i+m];
  256.       }
  257.       for(i=uend-m;i<=uend;i++) {
  258.          spc[i] = 0.0;
  259.          err[i] = 0.0;
  260.       }
  261.       return(0);
  262.    }
  263.    if(checkopt(argc,argv,"-slt",z)) {    /* shift left with time */
  264.       m=atoi(z);
  265.       for(i=ubeg;i<=uend;i++) {
  266.          spc[i] = spc[i+m];
  267.          err[i] = err[i+m];
  268.          tim[i] = tim[i+m];
  269.       }
  270.       for(i=uend;i<=uend;i++) {
  271.          spc[i] = 0.0;
  272.          err[i] = 0.0;
  273.          tim[i] = 0.0;
  274.       }
  275.       return(0);
  276.    }
  277.    if(checkopt(argc,argv,"-com",z)) {    /* compress spectra */
  278.       m=atoi(z);
  279.       i=ubeg; j=ubeg;
  280.       while(i <= (umax-m)) {
  281.          x=0.0;
  282.          for(k=0;k<m;k++) x=x+spc[i++];
  283.          spc[j] = x/m;
  284.          err[j] = err[i] / sqrt((double)m);
  285.          tim[j] = tim[j] * m;
  286.          j = j + 1;
  287.       }
  288.       umax=j-1;
  289.       return(0);
  290.    }
  291.  
  292.    if(checkopt(argc,argv,"-mir",z)) {    /* mirror spectra */
  293.       m=atoi(z);
  294.       i=ubeg; j=umax;
  295.       while(i < j) {
  296.          x=spc[i];
  297.          spc[i]=spc[j];
  298.          spc[j]=x;
  299.          x=err[i];
  300.          err[i]=err[j];
  301.          err[j]=x;
  302.          j=j-1; i=i+1;
  303.       }
  304.       return(0);
  305.    }
  306.  
  307.    if(checkopt(argc,argv,"-bel",z)) {    /* cut spectrum smoothly*/
  308.       m=atoi(z);
  309.       i=ubeg; j=uend;
  310.       while(i < j) {
  311.          x =  ((float) i) / ((float) m);
  312.          x=cos(x*1.57079);
  313.          if(i>=m) x=0.0;
  314.          spc[i] = x * spc[i];
  315.          err[i] = x * err[i];
  316.          i=i+1;
  317.       }
  318.       return(0);
  319.    }
  320.          
  321. /* ----------------------- monaric functions ------------------ */
  322.  
  323.    if(checkopt(argc,argv,"-sin",z)) {
  324.       while(ubeg <= uend) {spc[ubeg] = sin(spc[ubeg]); ubeg++ ;}
  325.       return(0);
  326.    }
  327.    if(checkopt(argc,argv,"-cos",z)) {
  328.       while(ubeg <= uend) {spc[ubeg] = cos(spc[ubeg]); ubeg++ ;}
  329.       return(0);
  330.    }
  331.    if(checkopt(argc,argv,"-exp",z)) {
  332.       while(ubeg <= uend) {spc[ubeg] = exp(spc[ubeg]); ubeg++ ;}
  333.       return(0);
  334.    }
  335.    if(checkopt(argc,argv,"-ln",z)) {
  336.       while(ubeg <= uend) {spc[ubeg] = log(spc[ubeg]); ubeg++ ;}
  337.       return(0);
  338.    }
  339.    if(checkopt(argc,argv,"-sqr",z)) {
  340.       while(ubeg <= uend) {spc[ubeg] = spc[ubeg] * spc[ubeg]; ubeg++ ;}
  341.       return(0);
  342.    }
  343.    if(checkopt(argc,argv,"-sqrt",z)) {
  344.       while(ubeg <= uend) {spc[ubeg] = sqrt(spc[ubeg]); ubeg++ ;}
  345.       return(0);
  346.    }
  347.    if(checkopt(argc,argv,"-abs",z)) {
  348.       while(ubeg <= uend) {spc[ubeg] = fabs(spc[ubeg]); ubeg++ ;}
  349.       return(0);
  350.    }
  351.    if(checkopt(argc,argv,"-int",z)) {
  352.       while(ubeg <= uend) {spc[ubeg] = floor(spc[ubeg]); ubeg++ ;}
  353.       return(0);
  354.    }
  355.       
  356. }
  357.  
  358. readthis(s,spc,err,tim)
  359. char s[];
  360. float spc[],err[],tim[];
  361. {
  362. int   max;
  363.  
  364.    max=readspec(s,spc,err,tim,s);
  365.    return(max);
  366. }
  367.  
  368. getnums(z,spc,err,tim)
  369. char z[];
  370. float spc[],err[],tim[];
  371. {
  372. int i,n,max;
  373.  
  374.    max=0;
  375.    if((z[0]=='-') || (z[0]<='9')) {
  376.       spc[max++]=atof(z);
  377.    } else {
  378.       max=readthis(z,spc,err,tim); /* possible spectra array */
  379.    }
  380.    return(max);
  381. }
  382.  
  383.